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EXTRACTION METHOD FOR STOKES FLOW WITH JUMPS IN 

THE PRESSURE 

KWANG SUNG CHANG AND DO YOUNG KWAK 

Abstract. In this paper, we consider a stationary, constant viscosity, in- 
compressible Stokes flow with singular forces along one or several interfaces. 
Assuming only the jumps of the pressure are present along the interface, we 
develop a new numerical scheme for such a problem. By constructing an ap- 
proximate singular function and removing it, we can apply a standard finite 
element method to solve it. A main advantage of our scheme is that one can 
use a uniform grid. We observe optimal O(h) order for the pressure and 0(h 2 ) 
order for the velocity. 

Key Words. Stokes equation, singular forces, jumps in the solution, extraction 
method, discontinuous pressure, uniform grid. 



1. Introduction 

In recent years, interface problems have become the subject of extensive re- 
search [3ll[T2[ll[T5l[T2l[ni[24l[26]. Many interesting physical phenomenons 
are described by the underlying partial differential equations having interface. For 
example, when two or more distinct materials or fluid with different conductivities, 
densities or permeability are involved, model equations often involves discontinuous 
coefficients to reflect the physical properties [TJ [5J [TU1 HU HH EH] • Often, the solu- 
tions of these interface problems must satisfy certain interface jump conditions due 
to physical conservation laws. Many numerical methods to deal with such prob- 
lems have been proposed [HI EU [22l [23] . Because of discontinuity of the solutions, 
standard numerical method do not yield accurate solutions even when fitted grid 
are used O [22l El] . 

In this paper, we propose an accurate, fast finite element method for Stokes 
problems with pressure jump conditions across a given interface which divides the 
domain into two parts. Any regular finite element meshes including uniform meshes 
are allowed. The idea is to consider certain singular function in a neighborhood of 
the interface whose jumps match the given jump conditions. By subtracting this 
function from the variational form, we obtain a new variational problem in which 
the solution has no jumps. To develop a numerical scheme, we choose one such sin- 
gular function and construct its approximation. The process consists of two main 
steps: First, we construct a piecewise linear function satisfying the jump conditions 
in a small strip near the interface. Then we extend it into whole region in some 
reasonable way. One of the natural method is to solve a harmonic/biharmonic equa- 
tion in one of the subdomains. It is quite natural to use finite element methods; 
However, the usual piecewise linear continuous finite element cannot be used be- 
cause of discontinuity of the data near the interface. Instead, Crouzeix-Raviart(CR) 
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Figure 1 . Sketch of the domain fl for the interface problem 



nonconforming element [31 H] which uses the midpoint of each edge as degree of 
freedom is appropriate in this case (see Section 4.1 for details). 

The next step is to subtract it from the variational formulation of Stokes problem 
which leads to standard variational form of Stokes problem. Some advantages of 
our scheme are: 

• We can use uniform mesh which is very efficient for moving interface such 
as time dependent problem. 

• Neither do we need adaptive mesh nor do we need extra degrees of freedom 
such as XFEM [22], yet our method achieve 0(h 2 ) for velocity and 0(h) 
for pressure, which is optimal with lowest order finite element. 

• The cost of constructing the discrete singular part is cheap since the work 
is equivalent to solving a Laplace problem in a subdomain. 

• After subtracting the singular part, the remaining task is equivalent to 
solving standard Stokes equation, hence it can be incorporated into the 
existing software. 

The outline of the paper is as follows: We start the formulation of the problem 
in Section 2, continued by a brief discussion of the removing singularity in the weak 
form of the Stokes equation. Section 3,4 describe the extraction method and its 
numerical scheme. In Section 5, some numerical results are shown. Conclusions 
follow in Section 6. 

2. Model Stokes problem 

Let O C M 2 be a convex polygonal domain. The domain fl is separated into 
two subdomains f2~ and fl + with £1 = f2~ U Cl + and f2~ n £1 + = 0. We assume 
that f2~ and f2 + are connected and <9f2~ PI <9f2 + = 0. The interface is denoted by 
T = Q~ nCl+. Let 



Then, we want to find the solution (u,p) e Hq (fi) 2 x Lq(Q) of the stationary homo- 
geneous Stokes problem for an incompressible viscous fluid confined in fl satisfies: 



(1) 



T 2 



{v e H 1 ^) : v = on dQ}, 
{qeL\n): f n q = 0}. 



(2a) 
(2b) 
(2c) 



/jAu + Wp = g + F, 
Vu = 0, 
u = 0, 



in f2, 
in f2, 
on d£l 
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with p = p(x) the pressure, u = u(x) the velocity, \i the constant viscosity, g G 
L 2 (fl) and F G i? _1 (fi), the external singular force. By [5],[B], this problem has 
a unique solution (u,p) G Hq(£1) 2 x Lq(Q,). The external singular force can be 
written as 



(F 1 (x,y),F 2 (x,y)) 



J /( S )£(x-X( S ))ds, 



where X(s) denotes the interface parameterized by s, / (s) is the force strength at 
this point, and S is the two-dimensional delta function. In general, this singular 
force leads to the jumps of the pressure and the velocity. However, dealing with 
those jumps for both the velocity and pressure is a heavy task, no one seems to 
have resolved it completely yet. Hence, in this paper, we restrict out attention to a 
rather simple case where the jumps are restricted the pressure only. So, we assume 
that the jumps of the pressure on the interface T are given by 

dp 



(3) 



Mr = J i( x )> 



<9n 



J 2 (x). 



Now, we define the related (affine) spaces. Let D + 
and define 

Hp(D) := H m (D+) n H m (D-), 
pmi(7i,T0(D) := jp g H^(D) [p] = 71, 



Dnn+, D~ 



with a piecewise norm 

Ml' 



dp 
dn 



= 72 on d n r 



{p G P m ;(7i,72)(£>)| p = on dD}, 



■(D) 



\\p\\ 2 h .. 



Ibll^ 



(D-y 



Here, for any domain in D, H m (D) is the usual Sobolev space of order m. 
ering the jump conditions, we decompose p as 



Consid- 



o 



(4) p = p 

where p° G H 1 (Q) D (^) and p* G V 1 '^ 1 ^ (O). In other words, p is splitted into 
regular part po and singular part p* . Using Green's theorem on each subdomains 
Q ± , we get 

(5) / Vp-vdx=— / pV-vdx+ / pv-nds. 

Jn± Jn± Jdn± 

Multiply a test function v in ([2a]) and integrating by part in each subdomain f2 + 
and f2~, we obtain a weak formulation of our problem as follows: find (u,p°) G 
H^(fl) 2 x Ll(Q) such that 

(6a) a(u,v) + b{v,p°) = (g, v) - b(v,p*)- < J x , v n > r , 
(6b) 6(u,q) = 0, 

where 



y vGH W, 

v ?ei 2 (!l), 



a(u, v):= / /iVu : Vv dx, 6(v, q) := — / qV-vdx, 
(Vu : Vv) := tr(VuVv). 

The resulting equation is the variational form of a standard Stokes equation with 
modified right hand side. The problem now is how to find an approximation to p* 
and approximate variational form. We will answer this question in the next section. 
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(a) Subdomains f2* , S in f2 . (b) fi^, QJ*, F in triangular grids. 

Figure 2. Some subdomains in Extraction method 



3. Extraction Method 

Taking the divergence of the equation (|2ap in each subdomain, we obtain the 
Poisson equation for the pressure: 



(7) 



- V 



Vg. 



We may split this equation into two equations 



(8a) 
(8b) 



Ap° = G 2 in il, 
Ap* = d in 



p° e H\n) nig(n), 



l;(Jl,Ja) 



(n), 



for some Gi and G2. Since p* has jump across T, the equations (|8b|) hold in f2 + 
and respectively, while (f8a|) holds in il. But solving this system numerically 
is not an easy task, since such a splitting is not unique and furthermore, p* has 
two interface conditions. In this paper, we propose the following method: First, 
we consider a narrow strip contained in whose outer boundary coincides with 



r (see Figure 2(a) I. Choose any function satisfying jump conditions and restrict 
it to S^call it p^ . Then we extend it into the whole Or by solving the following 
equation 



(9) 



Ap* = in fT \S =: fT, 

p* = P E on an*. 



Finally, set p* = on Q + . We call this scheme an Extraction Method(EM). The 
remaining task is to find a finite dimensional approximation to p* . 



4. Numerical Scheme 

Let Th be the usual quasi-uniform finite element triangulations of the domain Q. 
For any element T G 7^ , we call an element T an interface element if the interface 
r passes through the interior of T, otherwise we call T a non-interface element 
and we call an edge e G dT an interface edge if the interface T passes through the 
interior of e, otherwise we call e a non-interface edge. Now, we introduce some 
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(a) The numerical interface f\ (b) A typical interface triangle T 



Figure 3. Interface elements in the triangular mesh. 



notations: 



d I T 


= the set of all interface edges of an element T, 


T 1 


= the set of all interface elements, 


n-N 
J h 


= the set of all non-interface elements, 


nj, 


= U T > 








- \Jt-, 




TeT h 


K 


= n;-\n][. 


Even though the interface T is a curve in general, we replace the part of interface in 


T by the line segment connecting the intersection points with dT . Therefore, the 


interface T is replaced by its polygonal approximation T^. Henceforth, T is always 



assumed to be r^(see Figure 3(a) ). 



4.1. Construction of p* h . We construct p* h in two steps. First, we will consider 
p* h in il^. Suppose T is an interface element. For simplicity, we assume the three 
vertices are given by A\ = (0, 1), A 2 = (0, 0), ^3 = (1,0) (see Figure [3(b)] ) . For any 
element T in general position, all the constructions to be presented below carries 
over through afhne equivalence. Assume the interface meets with the element's 
edges at points B\ and B 2 - 

Let £i be the usual linear Lagrange nodal basis function associated with the 
vertex A, for i = 1, 2, 3. Then l\ = y, £2 = 1 — x — y, £3 = x. In the domain Q^, 
we write tp* as the following form: 



(10) r 



ip* — ol\£\ + a 2 £2 + ^£3 in T 
%//*+ = in T H 



Now imposing the jump conditions on T, we have: 
(11a) r~(Bi)-r + {Bi) = Jx{B x ), 

(lib) 1 p*-(B 2 )-x/ J *+(B 2 ) = JiGBa), 

(no arz^j^ _ 

on on 
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where B Q is the midpoint of B\B 2 . Then we have three unknowns CKj, i = 1,2, 3 in 
three equations (fTTa .b.c). Thus we can find coefficients {0^=1,2.3}- Note that two 
end points of the line segment T seg are located on the interface T, and hence the 
interface condition [p](x) = </i(x) is enforced exactly at these two end points, i.e., 
the point jump conditions (jlla.b) gives 

(l-a)a 2 +aa 3 = Ji(-Bi), 

bai + (1 - b)a 2 = Ji(B 2 ). 

The second condition on the interface segment r seg is the flux continuity. Hence, 
the derivative jump condition (|1 lc[) becomes 

(an Wi + a 2 V£ 2 + a 3 V£ 3 ) ■ = J 2 (B ). 

Since these conditions are represented by the following matrix equation: 



(12) 



1 — a a 
b l-b 

a —a — b b 



ai\ ( Ji(Bi) 
a 2 = Ji(B 2 ) 




"3 / V ^2(^0) 

the coefficient a^s of t/f* arc determined by the following formula: 

f at = (b(b-l) n + (a 2 +b) 32 -a(-b-l)p)/K 

(13) r : J a 2 = (b 2 n + a(aj2-bp))/K 

{ a 3 = ((a + 6 2 )j 1 + (a-l)(aj 2 ~6p))/X 

where K = a 2 + b 2 . We do this for every T £ n[. Having constructed pf on 
£1^ (lightly shaded region in Figure |2(b)[ ) , we now need to extend it into £l* h (dark 
shaded region in Figure 2(b) ), which will be done by solving the Laplace equation 

(1 4 ) ., _ ..L 



Ap* = in O*, 

p* = vl ondn* h , 



numerically. Numerical methods for this problem are well known; However, due to 
the discontinuity in the boundary data p^ on dfl^ , the usual nodal finite element 
space cannot be used. Instead, one can use Crouzeix-Raviart Pi -nonconforming 
finite element space [5] where the linear basis function has the degrees of freedom 
at the midpoints of edges. We denote it by i\. The finite element solution of (fl4|) 
on fi£ is denoted by p* h . Together with the construction above on we have 
obtained p* h in f2^ . 

4.2. Variational form after removing p* h . In this section, we explain how to 
remove the discrete singular part p* h from the weak formulation ([6]). First of all, we 
define some discrete spaces: 

Ph -~ iP e ^2(^/i)| p\t is constant for every T € T^}, 
P[ := { P eL 2 p\ T ± is linear for every T e T{ } , 

P h (D) := {p e L 2 (D) \ p is linear for every TnD, T G T h , 

p is continuous at the midpoints of the triangle edges} , 

V h := (P h (n h )) 2 . 

We now replace p* in ([6]) by p* hl but not without caution: Since p* h is now discon- 
tinuous along edges of T, we include line integrals. Thus we replace ((6aj) by the 
following form 

(15) a(u,v) + &(v,p°) = (g,v)-b(v,p* h ) + J(v,p* h ), 
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where 

J(v,P* h ) ■=- < K],v-n> r + V / p* h v-nds. 

The resulting equation can be solved again by a standard finite element method 
for Stokes problem: The simplest and most natural finite element method in this 
setting is the Crouzeix-Raviart finite element pair V \ x Thus, we have the 
following discrete Stokes problem: Find (u/j,p°) e Vh x P® such that 

(16a) a h (u h ,v h ) + b h (v h ,p° h ) = (9, v h ) - b h (v h ,p* h ) + J(v h ,p* h ), v v h eV h , 
(16b) b h (u h ,q h ) = 0, Veif, 

where 

a h (u h ,v h ) -=^2 ^ Vu/l : Vv '' dx ' b h (v h ,q h ) ■= ^V-v^dx. 



5. Numerical experiments 

In all of the experiments, the domain is a square and triangularized by uniform 
triangle grids with h x = h y = l/2 n ~ 1 for n = 3, • ■ • , 8. In order to describe 
the interface, we consider the level-set function $(x) for the interface T which is 
assumed to be smooth. Let $ : ft — > R be a continuous function such that 







' <0 


x in f2 _ 


(17) 


$(x) = < 


= 


x on r, 






> 


x in £1+ 



We assume that 3>(x) is smooth and V$ is not zero in any neighborhood of the 
interface T. Then the unit normal vector n(x) is represented by j^§r- The experi- 
ments in this subsection show that the method is robust. 

Example 1 (Constant Jump). The level-set function 3>(x), the jumps of the pres- 
sure and the boundary condition of the velocity are given as follows: 

$(x) = y/(x - 0.5) 2 + (y- 0.5) 2 - 0.25, 

J 150(a:- l/2)(y- 1/2), in VL+ 

\ 150(a:-l/2)(2/-l/2) + 30, mO- 
f Ul = -256a; 2 (a: - l) 2 y(y - l)(2y - 1), 
\ «2 = 256y 2 (y - l) 2 a;(a; - l)(2a; - 1), 
u = o?i dSl. 

with the domain Q = [0, 1] x [0, 1]. 



We observe the robust first order for the pressure and second order convergence 
for the velocity with L 2 -norm. 
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Table 1. (Constant jump) 1st order for the pressure and 2nd 
order for the velocity with || ■ \\ L2 



N x X Ny 


lb- Phil z,2 


Order 


II" - U h\\ L 2 


Order 


8x8 


5.1852 x 10-° 




1.6315 x 10" 1 




16 x 16 


2.0473 x 10"° 


1.34 


4.7486 x 10" 2 


1.78 


32 x 32 


8.1309 x 10- 1 


1.33 


1.2568 x 10~ 2 


1.92 


64 x 64 


3.6409 x lO" 1 


1.16 


3.2055 x 10~ 3 


1.97 


128 x 128 


1.7276 x lO" 1 


1.07 


8.0694 x lO" 4 


1.99 


256 x 256 


8.6181 x 10~ 2 


1.00 


2.0220 x 10~ 4 


2.00 




(a) p? in Q 




(b) p* in n{ 





(c) p* h in n (d) p h = p° h + p* in H 

Figure 4. pi, p^ in the constant jump case. 

Example 2 (Noncontant Jump). The level-set function $(x), the jumps 
pressure and the boundary condition of the velocity are given as follows: 

$(x) = y/x 2 + y 2 ~ 0.5, 

[p] r — 20 sin(x 2 y) — x 2 y, 



of 



dp 



On 



3x 2 y(20cos(x 2 y) - 1), 

!20xcos(a; 2 y)(2?/,x), in Q~ 
x(2y, x), in fi + 

on dft. 
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with the domain Q = [—1,1] X [—1,1]. 

We again have similar optimal convergence behavior. 

Table 2. (Nononstant jump) 1st order for the pressure and 2nd 
order for the velocity with || • || L2 



N x X Ny 


Ib-PhlL* 


Order 




Order 


8x8 


2.4718 x lCT 1 




1.6473 x 10" 2 




16 x 16 


9.5051 x 10~ 2 


1.38 


4.8611 x 10" 3 


1.76 


32 x 32 


4.2345 x 10" 2 


1.17 


1.5310 x 10~ 3 


1.67 


64 x 64 


2.0367 x 10^2 


1.06 


4.2185 x 10~ 4 


1.86 


128 x 128 


1.0152 x 10~ 2 


1.00 


1.0879 x 10~ 4 


1.96 


256 x 256 


5.0738 x 10~ 3 


1.00 


2.7698 x 10~ 5 


1.97 




6. Conclusions 

In this paper, we have introduced a new numerical method of solving Stokes 
interface problems having jumps in the pressure. The first step is to construct 
a piecewise linear function having small support in near the interface which 
satisfy the jump conditions. The second step is to extend it into by solving a 
discrete Laplace equation with Pi -nonconforming finite clement. Then removing it 
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from the original variational form, we obtain a Stokes problem with no jumps. The 
equation is then solved with the Crouzeix-Raviart nonconforming finite element 
pair. Our scheme is very effective since we can use any shape regular grid, not 
necessarily fitted grid. We have provided some numerical examples which show the 
optimal 0(h 2 ) error for velocity and 0(h) for pressure. 
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1. First section 

This is the first section. 

2. Second section 
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2.2. The second subsection. This is the first subsection of the second section. 
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ond subsection. 
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